Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix floating point error in bezier bounding box calculation #219

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Fix floating point error in bezier bounding box calculation #219

wants to merge 2 commits into from

Conversation

bouhek
Copy link

@bouhek bouhek commented Feb 1, 2024

Description

During my recent work with this library, I encountered an issue related to the calculation of bounding boxes for certain cubic Bezier curves. This problem appeared to stem from a floating point error that occurs during the computation of the denominator, which then affects subsequent conditions.

To resolve this, I have implemented a solution in this PR: the introduction of a FLOAT_EPSILON constant. This constant serves as a threshold for acceptable error in floating point calculations, effectively addressing the issue at hand.

To demonstrate the efficacy of this fix, I've added the zero_denominator_bezier_curve test. Included below are before-and-after images to illustrate the improvement. In these images, the Bezier curve is depicted by a black line, while the red rectangle represents the bounding box as calculated. This visual comparison clearly shows the enhancement brought about by this PR.

Before the fix

image

After the fix

image

Other test cases added

Test start_end_bbox_bezier_curve

image

Test general_bezier_curve

image

@bouhek
Copy link
Author

bouhek commented Feb 5, 2024

Hey, @mathandy, could you have a look at this when you have time? Thanks!

@Brishen
Copy link

Brishen commented Feb 13, 2024 •
edited
Loading

I found this error as well.
I would propose this method instead to avoid setting a constant.

from math import ulp

def bezier_real_minmax(p):
    """returns the minimum and maximum for any real cubic bezier"""
    local_extremizers = [0, 1]
    if len(p) == 4:  # cubic case
        a = [p.real for p in p]
        denom = a[0] - 3*a[1] + 3*a[2] - a[3]
        if abs(denom) < ulp(a[0]) and abs(denom) > 0:
            delta = a[1]**2 - (a[0] + a[1])*a[2] + a[2]**2 + (a[0] - a[1])*a[3]
            if delta >= 0:  # otherwise no local extrema
                sqdelta = sqrt(delta)
                tau = a[0] - 2*a[1] + a[2]
                r1 = (tau + sqdelta)/denom
                r2 = (tau - sqdelta)/denom
                if 0 < r1 < 1:
                    local_extremizers.append(r1)
                if 0 < r2 < 1:
                    local_extremizers.append(r2)
            local_extrema = [bezier_point(a, t) for t in local_extremizers]
            return min(local_extrema), max(local_extrema)

    # find reverse standard coefficients of the derivative
    dcoeffs = bezier2polynomial(a, return_poly1d=True).deriv().coeffs

    # find real roots, r, such that 0 <= r <= 1
    local_extremizers += polyroots01(dcoeffs)
    local_extrema = [bezier_point(a, t) for t in local_extremizers]
    return min(local_extrema), max(local_extrema)

@bouhek
Copy link
Author

bouhek commented Feb 13, 2024 •
edited
Loading

I found this error as well. I would propose this method instead to avoid setting a constant.

from math import ulp

def bezier_real_minmax(p):
    """returns the minimum and maximum for any real cubic bezier"""
    local_extremizers = [0, 1]
    if len(p) == 4:  # cubic case
        a = [p.real for p in p]
        denom = a[0] - 3*a[1] + 3*a[2] - a[3]
        if abs(denom) < ulp(a[0]) and abs(denom) > 0:
            delta = a[1]**2 - (a[0] + a[1])*a[2] + a[2]**2 + (a[0] - a[1])*a[3]
            if delta >= 0:  # otherwise no local extrema
                sqdelta = sqrt(delta)
                tau = a[0] - 2*a[1] + a[2]
                r1 = (tau + sqdelta)/denom
                r2 = (tau - sqdelta)/denom
                if 0 < r1 < 1:
                    local_extremizers.append(r1)
                if 0 < r2 < 1:
                    local_extremizers.append(r2)
            local_extrema = [bezier_point(a, t) for t in local_extremizers]
            return min(local_extrema), max(local_extrema)

    # find reverse standard coefficients of the derivative
    dcoeffs = bezier2polynomial(a, return_poly1d=True).deriv().coeffs

    # find real roots, r, such that 0 <= r <= 1
    local_extremizers += polyroots01(dcoeffs)
    local_extrema = [bezier_point(a, t) for t in local_extremizers]
    return min(local_extrema), max(local_extrema)

Hey @Brishen, thanks for the suggestion!
I like the idea of getting rid of the constant, but I feel like your solution wouldn't work for the arguments where denominator actually equals to 0 without the floating point error due to the and abs(denom) > 0 in your condition is bigger than ulp(a[0]). I'll try to revisit this tomorrow and add a test case for this too.

@Brishen
Copy link

Brishen commented Feb 13, 2024

@bouhek I think I got my conditions slightly wrong. It should probably be abs(denom) > ulp(a[0]) (and a[x] should probably be chosen to be the largest value).

from math import ulp

def bezier_real_minmax(p):
    """returns the minimum and maximum for any real cubic bezier"""
    local_extremizers = [0, 1]
    if len(p) == 4:  # cubic case
        a = [p.real for p in p]
        denom = a[0] - 3*a[1] + 3*a[2] - a[3]
        if abs(denom) > ulp(a[0]):
            delta = a[1]**2 - (a[0] + a[1])*a[2] + a[2]**2 + (a[0] - a[1])*a[3]
            if delta >= 0:  # otherwise no local extrema
                sqdelta = sqrt(delta)
                tau = a[0] - 2*a[1] + a[2]
                r1 = (tau + sqdelta)/denom
                r2 = (tau - sqdelta)/denom
                if 0 < r1 < 1:
                    local_extremizers.append(r1)
                if 0 < r2 < 1:
                    local_extremizers.append(r2)
            local_extrema = [bezier_point(a, t) for t in local_extremizers]
            return min(local_extrema), max(local_extrema)

    # find reverse standard coefficients of the derivative
    dcoeffs = bezier2polynomial(a, return_poly1d=True).deriv().coeffs

    # find real roots, r, such that 0 <= r <= 1
    local_extremizers += polyroots01(dcoeffs)
    local_extrema = [bezier_point(a, t) for t in local_extremizers]
    return min(local_extrema), max(local_extrema)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers
No reviews
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.
2 participants
@bouhek @Brishen